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ABSTRACT 

Rarefied  gas  phenomena  are  found  in  a  wide  variety  of  hypersonic  flow  applications,  and  accurate 
characterization  of  such  phenomena  is  often  desired  for  engineering  design  and  analysis.  In  atmospheric 
flows  involving  high  Mach  numbers  and/or  low  densities,  underlying  assumptions  in  a  continuum  fluid 
dynamics  description  tend  to  break  down,  and  distributions  of  both  molecule  velocity  and  internal  energy 
modes  can  depart  significantly  from  equilibrium.  In  this  paper,  characteristics  of  rarefied  and 
nonequilibrium  gas  flows  are  discussed,  and  several  numerical  methods  for  rarefied  flow  simulation  are 
outlined.  A  variety  of  hypersonic  flow  applications  are  presented,  and  additional  discussion  is  provided 
for  directions  of  current  and  future  research. 


1.0  INTRODUCTION 

The  field  of  rarefied  gas  dynamics  is  concerned  with  gas  flow  problems  that  do  not  conform  to  a 
continuum  description;  in  particular,  flows  for  which  the  Knudsen  number  -  defined  as  the  ratio  of  the 
mean  free  path  to  some  characteristic  length  -  is  greater  than  that  required  to  approximate  the  flow  as  a 
continuous  medium.  The  mean  free  path  is  in  turn  defined  as  the  average  distance,  typically  in  a  reference 
frame  which  moves  with  the  flow,  that  a  molecule  travels  before  colliding  with  another  molecule. 
Traditional  continuum  gas  dynamics  models  tend  to  break  down  due  to  gradient  transport  assumptions  for 
diffusion  terms  in  continuum  equations  of  mass,  momentum,  and  energy  conservation.  More  generally,  in 
a  rarefied  flow,  the  effect  of  intermolecular  collisions  to  drive  the  distribution  of  molecular  velocities 
toward  the  equilibrium  limit  (called  a  Maxwellian  velocity  distribution)  is  insufficient  to  maintain  a  near¬ 
equilibrium  distribution.  As  a  result,  equilibrium  assumptions  underlying  the  continuum  Navier-Stokes 
gas  dynamics  conservation  equations  become  invalid,  and  detailed  characterization  of  the  velocity 
distribution  function  is  required  for  accurate  numerical  simulation  of  the  flow. 

The  above  characteristics  are  typically  found  where  length  scales  of  gradients  for  the  macroscopic 
quantities  are  on  the  order  of  the  molecular  mean  free  path  of  the  gas.  If  the  Knudsen  number  is  used  to 
characterize  the  degree  of  rarefaction  effects,  then  we  can  expect  rarefied  flow  to  result  from  either  a  large 
mean  free  path  or  small  characteristic  length  scales.  Gas  flow  problems  involving  large  mean  free  paths 
are  often  encountered  in  high  altitude  atmospheric  flows,  such  as  those  around  reentry  vehicles  and 
spacecraft  in  low  Earth  orbit;  vacuum  systems  as  used  in  industrial  processes  including  chemical  vapor 
deposition  and  freeze  drying;  and  in  low  density  plume  flows  from  spacecraft  or  rockets.  Flow  problems 
for  which  rarefaction  effects  occur  due  to  small  length  scales  include  hypersonic  atmospheric  flows,  where 
local  characteristic  flow  lengths  may  be  very  small  inside  thin  shock  and  boundary  layer  regions;  low 
speed  flows  through  or  around  micro-electro-mechanical  systems  (MEMS);  and  other  mechanical  devices 
-  such  as  computer  hard  drives  -  which  involve  small  air  gaps. 

Early  applications  of  rarefied  gas  dynamics  involved  problems  in  vacuum  technology,  but  the  field  has 
had  a  major  impact  on  aerothermal  analysis  and  design  of  spacecrafts,  satellites,  missiles,  etc.  flying  at 
high  altitudes  [1,  2].  In  these  applications,  rarefaction  effects  become  important  due  to  very  low  densities 
and/or  very  high  Mach  numbers.  In  recent  years,  however,  several  applications  have  emerged  for  which 
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ORGANIZATION 


strong  rarefaction  effects  occur  as  a  result  of  small  length  scales  associated  with  the  flowfield  geometry; 
one  example  of  this  is  the  MEMS  device.  In  this  paper,  however,  we  will  be  primarily  concerned  with  the 
applications  of  rarefied  gases  in  hypersonic  flows. 


2.0  RAREFIED  FLOW  PHENOMENA 

2.1  Mean  Free  Path  and  Various  Flow  Regimes 

Rarefied  gas  flows  are  primarily  classified  in  three  flow  regimes: 

•  0.01  <  Kn  <  0.1  Slip  flow  regime 

•  0.1  <  Kn  <  10  Transitional  regime 

•  Kn  >  10  Free  molecular  regime 

Applications  in  slip  flow  regime  manifest  mainly  near  the  boundaries  of  vehicles  in  high  speed  flight,  as 
the  so  called  velocity  slip  and  temperature  jump  phenomena,  which  are  treated  by  special  boundary 
conditions.  For  the  rest  of  the  flowfield,  the  continuum  gas  dynamics  equations  are  still  valid.  For  the 
other  extreme  of  the  free  molecular  regime  where  Kn  is  high,  the  gas  is  very  rarefied.  In  this  regime, 
where  the  mean  free  path  is  much  larger  than  the  characteristic  dimension  of  the  body,  the  collisions  occur 
infrequently  even  after  the  particles  hitting  the  surface  reflect  away.  In  the  transitional  regime  (0.1  <  Kn  < 
10),  both  the  collisions  between  particles  and  with  the  surface  assume  importance,  hence  this  is  an  area  of 
rarefied  gas  dynamics  with  considerable  body  of  research.  The  primary  thrust  of  this  paper  is  for  flows  in 
the  transitional  regime. 

In  the  molecular  description  of  a  gas,  the  kinetic  theory  is  based  on  the  hypothesis  that  the  gas  is 
composed  of  molecules.  Although  for  a  gas  at  rest,  the  molecules  are  separated  by  distances  large  enough 
that  they  interact  weakly  with  each  other.  But  when  the  molecules  interact  strongly  they  are  considered  to 
be  collisions.  In  the  treatment  of  interactions,  if  the  time  average  of  potential  energy  can  be  neglected 
compared  with  the  kinetic  energy,  the  gas  is  called  perfect.  When  this  condition  is  not  satisfied,  the  gas  is 
called  van  der  Waals  gas.  The  vast  body  of  the  literature  on  rarefied  gas  dynamics  deals  with  gases  that  are 
perfect.  In  practice,  gases  of  neutral  molecules  up  to  pressures  of  hundreds  of  atmospheres  are  considered 
perfect.  Another  useful  way  to  classify  the  regime  is  with  pressure.  For  gases  even  up  to  very  high 
pressures,  three-body  collisions  are  small  in  comparison  with  two-body  collisions,  and  here,  Classical 
Newtonian  mechanics  is  used  to  describe  the  motion  of  the  molecules  except  where  quantum  effects 
become  important.  In  addition  to  the  importance  of  quantum  effects  at  very  low  temperatures  and  for  light 
molecules,  which  we  will  not  consider,  quantum  effects  must  be  taken  into  account  in  aerospace 
applications  for  inelastic  collisions  of  atoms  and  molecules  where  the  internal  degrees  of  freedom  are 
excited.  The  potential  of  these  molecular  interactions  may  be  calculated  only  be  the  use  of  quantum 
mechanics. 

2.2  Breakdown  of  Navier-Stokes  Equations 

The  Navier  Stokes  equation  are  derived  assuming  only  small  deviations  from  equilibrium  as  characterized 
by  a  Maxwellian  distribution  function.  In  the  equilibrium  state  of  a  gas,  the  distribution  of  molecular 
energy  states  does  not  vary  with  time.  Although  intermolecular  collisions  occur,  and  energy  is  exchanged 
in  collisions,  the  gas  maintains  the  same  number  of  particles  in  a  given  energy  state.  Similarly,  if  the  gas 
were  to  undergo  a  reaction,  the  equilibrium  state  would  imply  that  the  rate  at  which  the  forward  reaction 
progresses  is  exactly  equal  to  the  rate  at  which  the  reverse  reaction  progresses,  such  that  the  composition 
of  the  gas  does  not  vary  over  time.  On  the  other  hand,  a  state  of  nonequilibrium  is  often  a  transient  state;  it 
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is  the  process  by  which  the  gas  approaches  a  new  equilibrium  state. 

The  continuum  equations  are  also  invalid  in  regions  of  very  low  density.  When  a  relatively  small  number 
of  particles  occupy  the  flowfield  of  interest,  the  continuum  hypothesis  is  obviously  invalidated  as  large 
regions  of  the  flowfield  may  contain  no  particles  at  all.  This  may  occur  naturally  as  in  the  upper 
atmosphere,  or  may  be  caused  by  rapid  expansion  of  the  flow.  In  this  case,  the  collision  rate  drops 
drastically  and  the  continual  decrease  in  temperature  predicted  by  the  continuum  equations  cannot  be 
maintained.  This  is  because  as  the  number  of  collisions  experienced  by  the  particlescontinues  to  decrease, 
the  translational  temperature  eventually  levels  out  or  “freezes”  as  the  free  molecular  regime  is  approached 
[3,  4,  5].  Alternatively,  when  the  length  scale  of  consideration  is  on  the  order  of  the  mean  free  path  in  the 
gas,  the  flow  may  “appear”  rarefied  even  at  normal  densities,  and  the  continuum  hypothesis  will  no  longer 
hold.  This  case  is  of  practical  interest  for  the  analysis  of  microscale  aerodynamics  and  fluid  flow  through 
micro-electro-mechanical  systems  (MEMS)  [6]. 

The  validity  of  the  continuum  equations  is  inherently  tied  to  the  collision  rate  in  the  gas.  Letting  v 
represent  the  intermolecular  collision  rate  and  tf  represent  the  characteristic  flow  time,  the  Knudsen 
number  Kn  may  be  defined  as  follows. 


Kn  =  —  (1) 


Alternatively,  letting  c  be  the  mean  molecular  speed,  X  the  mean  free  path  (the  average  distance  traveled 
by  a  molecule  before  encountering  a  collision),  and  L  represent  a  characteristic  length  scale,  the  above 
equation  becomes 


vtf 


1 

(< c/X)(L/c ) 


A 

L 


(2) 


The  validity  criterion  for  the  continuum  equations  has  been  stated  as  Kn  «  1.  However,  the  choice  of 
characteristic  length  scale  may  be  somewhat  ambiguous.  Choosing  a  length  scale  based  upon  the  geometry 
over  which  the  flow  is  considered  results  in  a  parameter  which  may  describe  globally  how  well  the 
continuum  equations  apply,  however,  it  does  not  address  the  properties  of  the  local  flow  physics  which 
may  occur.  A  better  choice  is  to  consider  a  local  length  scale  related  to  the  physics  occurring  in  the 
flowfield.  This  may  be  a  quantity  such  as  a  shock  thickness,  boundary  layer  thickness,  or  some  other 
length  associated  with  the  local  flow  physics. 

The  limits  shown  in  Figure  1  are  based  upon  experimental  observations  and  should  not  be  thought  of  as 
hard  limits  or  having  a  rigorous  mathematical  basis.  Indeed  this  is  the  very  reason  it  is  desired  to  study  the 
criteria  for  continuum  breakdown.  The  Euler  equations  are  strictly  based  upon  the  Maxwellian  velocity 
distribution,  which  is  the  prevailing  distribution  when  the  gas  is  in  equilibrium.  Furthermore,  as  they 
include  no  viscous  terms,  they  are  only  valid  at  very  low  Knudsen  number,  or  when  the  collision  rate 
becomes  quite  large.  When  this  occurs,  the  gas  is  able  to  redistribute  its  energy  to  accommodate  varying 
conditions  almost  instantly,  and  the  flow  is  practically  in  equilibrium  everywhere. 

The  Navier-Stokes  equations  exhibit  an  extended  Knudsen  number  validity  over  the  Euler  equations  as 
seen  in  Figure  1 .  Pinpointing  exactly  where  they  become  invalid  is  somewhat  difficult  to  predict,  but  it  is 
generally  accepted  that  this  occurs  at  Knudsen  numbers  somewhere  around  0.1  [4].  The  Navier-Stokes 
equations  can  be  derived  by  using  the  first  Chapman-Enskog  correction  to  the  Maxwellian  distribution 
function  [7]  [8].  This  means  that  the  Navier-Stokes  equations  are  based  upon  a  small  departure  from  the 
equilibrium  distribution  function  and  will  fail  in  regions  of  significant  nonequilibrium. 
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Many  flow  scenarios  contain  several  physical  processes  that  exhibit  local  Knudsen  numbers  at  many 
different  locations  on  this  spectrum.  This  is  particularly  the  case  in  some  hypersonic  flows.  Consider  a 
typical  blunt  body  in  hypersonic  flow.  The  physics  involved  with  the  associated  structure  of  the  strong 
shocks  invalidate  the  continuum  equations.  Aft  of  the  shock  the  continuum  equations  may  hold.  Close  to 
the  body,  strong  gradients  exist  across  the  boundary  layer  and  the  process  of  activating  the  internal  modes 
of  the  gas  may  occur  in  nonequilibrium  fashion.  Further  downstream  the  flow  may  expand  around  the 
body  to  the  point  at  which  the  density  is  too  low  for  the  continuum  equations  to  hold.  Such  a  mixed 
flowfield  exemplifies  why  a  hybrid  method  may  be  desirable,  illustrates  the  importance  of  understanding 
where  the  continuum  equations  are  and  are  not  valid,  and  emphasizes  the  need  for  a  good  understanding  of 
nonequilibrium  phenomena. 

In  the  past,  several  parameters  have  been  examined  in  an  attempt  to  quantify  continuum  breakdown.  For 
obvious  reasons,  many  of  the  initial  investigations  considered  parameters  in  the  form  of  a  Knudsen 
number  based  upon  the  local  flow  physics.  This  approach  is  based  on  the  observation  that  if  flow  variables 
change  over  relatively  small  distances  in  comparison  to  the  local  mean  free  path,  it  is  unlikely  that  enough 
collisions  are  experienced  for  the  gas  to  reach  equilibrium  at  the  new  conditions.  Bird  [3]  was  one  of  the 
first  to  examine  the  validity  of  such  a  parameter  using  DSMC.  He  did  so  in  the  context  of  nonequilibrium 
resulting  from  expansions  of  the  gas.  In  doing  so,  he  suggested  the  following  parameter  as  a  means  of 
quantifying  equilibrium  breakdown. 


1 

D(lnp) 

V 

Dt 

(3) 


where  p  is  the  density  of  the  gas.  In  the  case  of  one  dimensional  steady  flow,  this  becomes 


P 


dp 

dp 

dx 

V  8  p 

dx 

(4) 


where  u  is  the  bulk  gas  velocity  and  M  is  the  Mach  number  in  the  gas.  Bird  found  that  a  value  of  P  greater 
than  0.02  tended  to  predict  the  onset  of  nonequilibrium  well. 

While  this  parameter  seems  to  indicate  continuum  breakdown  well  in  the  specific  case  of  gaseous 
expansions,  it  does  not  predict  nonequilibrium  well  in  regions  of  low  velocity.  As  seen  in  the  above 
equation,  the  Mach  number  dependence  will  drive  P  to  zero  in  regions  of  low  velocity  no  matter  what  the 
degree  of  nonequilibrium.  Boyd  proposed  examining  a  Knudsen  number  based  upon  the  gradient  local 
length  (GLL)  as  a  means  of  examining  continuum  breakdown  [9]. 

KnGLL=^\VQ\  (5) 


where  Q  is  some  flow  property  of  interest.  The  dependence  of  Bird's  parameter  on  a  gradient  based  local 
Knudsen  number  is  readily  seen  in  the  following  equation. 


P  =  M 


(6) 


where  Kn  = 


dp 


dx 


is  the  Knudsen  number  based  upon  the  gradient  local  length  of  the  density. 


Determining  which  flow  variable  is  the  best  to  use  for  examining  continuum  breakdown  is  nontrivial,  and 
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may  vary  from  one  flow  scenario  to  another.  For  this  reason,  Wang  and  Boyd  [10]  proposed  using  the 
following  criterion 


Knmia  =  max  [  Knp ,  Kitr ,  Kriy  }  (7) 

where  Knr  is  the  Knudsen  number  based  upon  the  gradient  local  length  of  the  temperature,  and  Knv  is  the 
Knudsen  number  based  upon  the  gradient  local  length  of  the  fluid  velocity.  They  proposed  that  significant 
nonequilibrium  was  indicated  when  this  parameter  became  greater  than  about  0.05 

Unfortunately,  while  such  parameters  are  relatively  easy  to  compute,  they  do  not  necessarily  have  a  firm 
theoretical  grounding.  Qualitatively,  they  make  some  sense,  in  that  many  forms  of  nonequilibrium  are  set 
into  motion  by  the  presence  of  strong  gradients  in  the  flow  variables.  Such  parameters  are  also  tied  to  the 
traditional  understanding  of  a  Knudsen  number  restriction  on  the  continuum  equations.  However,  the  fact 
that  multiple  parameters  must  be  computed  suggests  that  the  continuum  prediction  criterion  is  problem 
dependent  at  this  stage  of  development.  It  should  be  noted  that  the  computation  of  the  heat  flux  and  shear 
stress  terms  involves  the  combination  of  a  number  of  flow  gradients.  This  partially  speaks  to  why  multiple 
gradient-based  Knudsen  parameters  are  required  to  obtain  a  reasonable  predictive  capability. 

2.3  Nonequilibrium  Internal  Energy  Excitation 

Thermal  nonequilibrium  among  the  internal  modes  of  motion  of  the  molecular  species  in  supersonic  and 
hypersonic  flows  has  major  influences  on  flowfield  temperatures,  radiation  signatures,  and  aerodynamic 
heating.  Prediction  of  such  flowfields  by  computer  modeling  techniques  requires  a  detailed  knowledge  of 
the  rates  and  mechanism  of  energy  storage  and  transfer  among  the  internal  modes  of  molecular  motion. 
Among  the  large  body  of  literature,  we  mention  a  few:  Anderson  [11],  Capitelli  [12],  Chernyi  et  al.[13], 
Rich  et  al.  [14],  and  Adamovich  et  al.  [15]. 

The  various  modes  of  motion  for  diatomic  gas  molecules  which  are  of  interest  in  aerothermodynamics  are 
translation,  rotation,  vibration,  electronic  and  the  mutual  interaction  of  these  modes.  In  addition,  molecules 
can  dissociate  and  ionize.  For  our  purposes  of  engineering  interest,  translational  motion  can  be  described 
by  classical  mechanics  and  is  treated  in  other  sections  in  this  paper.  However,  to  quantify  the  degree  of 
nonequilibrium  in  the  internal  energy  modes  one  would  have  to  consider  the  internal  energy  modes  of 
rotation,  vibration,  and  electronic  excitation,  all  of  which  are  quantized.  When  the  gas  is  in  a  state  of 
thermal  equilibrium,  the  distribution  of  internal  energy  among  the  various  modes  is  said  to  be  governed  by 
a  Maxwell-Boltzmann  distribution,  which  can  be  written  as  the  fractional  number  of  molecules  rii  in  the  ith 
internal  energy  state  as: 


n  a 


-exp 


EtM 

kT 


where  the  internal  partition  function  Qint  is  given  by 


Qint  =E‘?-eXP 


C  E  ^ 

I~‘i,  int 

~w 


(8) 


(9) 


where  N  =  ^n,  is  the  total  number  of  atoms  and  molecules  and  g,  is  the  statistical  weight  of  the  ith 

i 

internal  energy  state,  k  is  the  Boltzmann  constant,  and  T  is  the  kinetic  temperature  of  the  gas.  In  this 
section,  we  will  review  some  recent  advances  concerning  situations  where  the  rate  of  relaxation  to  the 
equilibrium  distribution  given  by  Eq.  (8)  is  slow  relative  to  the  mean  collision  time.  Each  of  the  energy 
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modes  must  undergo  a  finite  number  of  molecular  collisions,  governed  by  a  relaxation  time,  in  order  to 
equilibrate  with  the  translational  mode.  The  relative  magnitude  of  the  relaxation  times  may  be  stated  as 
follows: 


T0~Tt<  Tr  «  TV  <  Td 


where  r0  is  the  mean  collision  time  and  tt,  rR ,  zy,  and  rd  are  respectively  the  relaxation  times  for 
translational,  rotational,  vibrational  energy  modes  and  dissociation. 

For  a  diatomic  molecule  N2 ,  the  symbolic  equations  governing  the  Vibrational-Translational  (VT) 
transitions  responsible  for  the  variation  of  the  particles  distributed  in  the  ith  vibrational  level  are 

N2(i)  +  N2^±N2(i')  +  N2  (10) 

and  the  equations  governing  the  Vibration- Vibration  (VV)  process  are 

N2  ( i )  +  N2  ( j )  ^  N2  (f )  +  N2(f)  (11) 

The  equations  governing  the  dissociation  process,  considering  dissociation  to  take  place  from  any 
vibrational  level,  are 


N2(i)  +  M  ^2N  +  M 


(12) 


The  kinetics  of  vibrational  energy  exchanges  among  the  quantum  states  are  simulated  using  the  vibrational 
master  equations  [16],  for  which  the  population  distributions  are  calculated  with 


COy  = 


X \}-vT  (z '  — >  i ) pv  p  ~  kyj  (z  — : >  i ') pv,p\ 

V 

+  z  Lrw  U  j) Pi  Pj'  -  rw  (i,  j  ->  i j ') PtPj ] 

+  \_~Kd  {i  continuum) ptpNi  ~kVD  (continuum  —>i)pNpN^ 


(13) 


where  is  the  molecular  weight  and  p  indicates  the  relative  population  at  a  given  state.  Energy 
exchanges  consist  of  multiquantum  transitions  for  VT  and  VV  processes.  Additionally  dissociation  can 
occur  from  any  of  the  vibrational  states  to  continuum.  The  VT  process  is  associated  with  the  rate 
coefficient  kyr  where  the  molecule  loses  or  gains  a  number  of  vibrational  quanta.  The  de-excitation  rate 
from  i  ’  to  i  for  colliding  molecules  is  denoted  by  kV7{i  the  inverse  collision  from  i  to  i  ’  by  kyiii^i  0- 
When  considering  VV  exchanges,  the  initial  and  final  vibrational  states  of  each  collision  partner  must  be 
identified;  thus  the  transition  rate  from  V  to  i  and  j’  to  j  is  given  by  rvv(i,J,—>i,j).  Consistency  of  the  rate 
coefficient  with  the  principle  of  detailed  balance  is  enforced. 

The  vibrational  energy  of  the  diatomic  molecules  /V2,  treated  as  anharmonic  oscillators,  is  given  in  terms 
of  the  quantum  level  energies  and  state  densities  by 


e  =y^£ 

cvib 

i=o  P 


(14) 


In  this  equation,  pjp  is  the  fractional  population  in  the  ith  vibrational  level  and  is  the  quantum  level 


1  -6 


RTO-EN-AVT-194 


Review  of  Rarefied  Gas  Effects  in  Hypersonic  Applications 


energy  given  by  the  third-order  approximating  formula 


he 


CQ„ 


V  ^ J 


(15) 


where  h  is  the  Planck’s  constant  and  c  is  the  speed  of  light.  Here,  sx  denotes  ground  state  vibrational 
energy,  s2  denotes  first  excited  state,  and  so  on.  For  nitrogen,  the  spectroscopic  constants  are  given  by 
coe  -  2358.57  cm'1,  coexe  =  14.324  cm'1,  and  coeye  =  -0.00226  cm'1.  When  /*  =  48,  the  vibrational  energy 
value  exceeds  the  dissociation  energy  of  9.86  eV. 

Vibrational  relaxation  involves  two  widely  different  time  scales,  associated  with  the  VT  and  VV 
processes.  The  relaxation  in  diatomic  gases  may  be  significantly  affected  by  the  rapid  VV  transfer 
processes.  In  VV  exchange  collisions  between  two  molecules,  the  vibrational  quanta  in  one  molecule 
increases  while  in  the  other  it  decreases.  Due  to  anharmonicity,  this  results  in  a  net  energy  transfer  with 
the  translational  mode.  Such  transfers  occur  with  high  probability  because  the  net  change  in  energy  of 
each  VV  transitions  is  small  relative  to  VT  transitions.  In  fact,  the  difference  in  time  scales  is  often 
tvv  «  Tyr  so  that  we  can  identify  the  VV  rates  as  “fast”  rates  and  the  VT  rates  as  the  “slow”  rates  in  a 
relaxation  process.  During  the  initial  stages  of  relaxation,  the  VV  exchanges  dominate  and  during  the 
remaining  bulk  of  relaxation,  the  VT  transitions  are  also  effective.  In  a  single  component  system  of 
anharmonic  oscillators,  VV  processes  may  have  significant  effects  for  all  quantum  levels  except  the  upper 
levels.  For  large  vibrational  quantum  numbers  and  closely  spaced  energy  levels  of  an  anharmonic 
oscillator,  the  VT  transitions  are  rapid.  The  population  distribution  of  these  upper  levels  is  normally  well 
described  by  a  vibrational  temperature  equal  to  the  translational  temperature. 

For  VT  processes  which  transfer  energy  between  the  vibrational  and  translational  degrees  of  freedom,  the 
energy  exchange  is  of  the  order  of  the  vibrational  energy  spacing.  VV  processes,  affected  by  the 
anharmonicity  of  the  system,  influence  the  flow  of  vibrational  energy  transfer  between  lower  and  upper 
levels.  We  define  Tvib  as  the  vibrational  temperature  based  on  relative  populations  at  all  vibrational  levels, 
and  Tfl  is  the  temperature  based  on  the  ratio  between  first  excited  state  and  ground  state  populations.  For 
Tfi  ^  Tyib  the  vibrational  distribution  differs  from  that  of  a  Boltzmann,  the  deviation  increasing  with 
increasing  vibrational  quantum  level.  For  Tfl  <  Tvib  anharmonicity  causes  the  underpopulation  of  upper 
vibrational  levels,  whereas  for  Tfl  >  Tvib  it  causes  an  overpopulation.  The  steady-state  nonequilibrium 
vibrational  distribution  function  in  a  one  component  system  of  anharmonic  oscillators  was  derived  by 
Treanor,  Rich  and  Rehm  [17].  Owing  to  the  high  vibration- vibration  energy  exchange  in  the  anharmonic 
oscillators,  the  low  to  intermediate  vibrational  levels  achieved  a  concave -upward  population  distribution 
known  as  the  Treanor  distribution.  A  popular  assumption  in  modern  computations  of  hypersonic 
aerodynamic  simulations  is  to  treat  the  molecule  as  a  simple  harmonic  oscillator  where  the  vibrational 
energy  is  modeled  by  the  Landau-Teller  method  [8]. 

The  rates  (cm3/s)  for  VV  and  VT  exchange  in  nitrogen  are  presented  in  Figure  2  at  6000K.  The  VT  rate 
corresponds  to  j+ 1  — »  j  increases  monotonically  with  the  quantum  level  vibrational  energy.  Two  relevant 
VV  exchanges  are  also  presented,  both  weighted  by  the  relative  equilibrium  population  of  the  vibrational 
state  of  the  collision  partner:  The  rate  of  weighted  VV  exchange  involving  the  ground  state  increases  with 
quantum  level  vibrational  energy  and  then  decreases  due  to  increasing  disparity  of  the  energy  being 
exchanged.  The  rate  of  population-weighted  resonant  exchange  increases  with  vibrational  quantum 
number  initially  and  then  decreases  due  to  the  reduced  population  associated  with  the  excited  state  of  the 
collision  partner. 

Also  shown  in  Figure  2  is  the  equilibrium  dissociation  rate.  This  is  the  value  required  to  support  the 
equilibrium  dissociation  rate  from  a  single,  given  vibrational  state.  In  the  "ladder  model"  the  state  of 
interest  is  the  last  bound  vibrational  state.  For  a  given  set  of  VV  and  VT  transition  rates,  the  relative 
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importance  of  the  VV  terms  r.^.  ,  the  VT  term  kj+ij  ,  and  the  dissociation  term  can  be  assessed.  It  can  be 

concluded  that  for  nitrogen  at  these  temperatures,  the  dissociation  probability  for  the  last  state  is  high 
relative  to  both  the  weighted  VV  and  VT  rates,  and  that  VT  is  more  important  than  VV.  Regarding 
depletion  of  the  dissociating  state,  one  would  anticipate  that  depletion  will  be  significant  for  Rate  Set  1 
(Figure  la)  [18,  19]  because  the  dissociation  rate  exceeds  the  V-T  rate  considerably  and  not  quite  as 
significant  for  Rate  Set  2  (Figure  lb)  [20].  Since  the  VT  probabilities  are  higher  than  those  of  weighted 
VV  energy  transfers,  the  effect  of  VV  on  the  population  depletion  is  negligible,  as  shown  in  Figure  2  [16]. 
The  aforementioned  case  is  a  good  example  of  the  need  to  have  reliable  transition  rates  for  accurate 
predictions  with  state  specific  models. 

Figure  3  illustrates  VT  multiquantum  and  single-quantum  rate  coefficients  from  the  25th  vibrational 
quantum  level  (i  =  25)  to  lower  levels  (/'  =  15,  23,  24)  [21,  22,  23].  The  multiquantum  rate  coefficients 
(/’  =15,  23)  become  comparable  to  that  of  single  quantum  (/'  =  24)  with  increasing  translational 
temperature.  It  may  be  concluded  that  multiquantum  processes  are  important  at  high  temperature 
conditions.  On  the  contrary,  large  multiquantum  transitions  are  less  probable  at  lower  temperatures  (see 
V  -  15  curve)  since  collision  energies  are  not  high  enough  to  permit  large  energy  jumps.  Because  the 
energy  levels  are  closer  to  each  other  with  increasing  vibrational  quantum  numbers,  the  vibrational  jumps 
are  more  probable  at  moderate  and  high  vibrational  levels.  Therefore,  one  may  expect  from  Figure  3  that 
the  lower  vibrational  levels  are  populated  mainly  by  single-quantum  transitions. 

Figure  4  shows  the  influence  of  multiquantum  VT  processes  on  the  evolution  of  the  translational 
temperature  behind  the  standing  shock  wave.  The  translational  temperature  due  to  single-quantum 
transitions  decays  much  more  gradually  with  distance.  One  can  see  that  at  the  higher  temperatures, 
immediately  behind  the  shock  front,  the  effect  of  multiquantum  transitions  is  very  high,  leading  to  a  faster 
rate  of  equilibration.  The  differences  in  macroscopic  quantities  from  single-  and  multiquantum  effects  can 
be  attributed  to  the  difference  in  the  total  transition  probability  of  a  VT  exchange;  i.e.  for  the  probability  to 
sum  up  to  unity,  all  multiquantum  transitions  have  to  occur.  This  constraint  would  imply  that  the 
assumption  of  single  quantum  energy  exchanges  will  result  in  a  total  probability  less  than  unity.  The 
difference  in  the  total  probability  of  energy  exchange  causes  the  mean  energy  transfer  to  be  different  for 
the  single  and  multiquantum  transitions,  as  seen  for  the  translational  temperature  difference  in  Figure  4. 

Figure  5  illustrates  the  vibrational  population  distribution  as  a  function  of  the  vibrational  energy  in  the 
quantum  energy  states,  for  multiquantum  and  single-quantum  processes  at  various  distances  (0.01  and 
0.03  cm)  behind  the  shock.  It  is  clear  that  multiquantum  mechanisms  populate  intermediate  and  high 
levels  more  quickly,  i.e.  for  i  >  10,  vibrational  energy  greater  than  3  eV.  For  the  multiquantum  transition 
case,  the  upper  vibrational  levels  follow  a  Boltzmann  distribution  at  the  higher  temperature.  However, 
concentrations  of  the  low  energy  levels  (i  <  10,  vibrational  energy  less  than  3  eV)  vary  and  decrease 
slowly  during  the  relaxation  of  the  gas  flow  behind  the  shock  wave.  The  slow  variation  occurs  because  VT 
transitions  in  the  lower  vibrational  states  are  required  to  overcome  a  large  energy  gap,  and  consequently 
single  quantum  and  multiquantum  collisions  are  less  effective.  The  concentrations  of  the  lower  vibrational 
levels  experience  a  lag  because  molecules  jump  to  higher  energy  states  due  to  VT  exchanges.  The  inset  in 
Figure  5  shows  that  the  population  is  smallest  in  the  low-lying  states  when  the  flow  reaches  equilibrium 
downstream  of  the  shock  wave.  In  the  vicinity  of  the  shock  wave,  at  the  locations  of  0.03  cm  and  0.01  cm, 
the  population  in  these  low-lying  states  is  higher  when  multiquantum  jumps  are  considered  and  highest 
with  single  quantum  jumps.  However,  with  an  increase  in  vibrational  levels,  the  equilibrium  population 
increases  to  a  higher  value  than  the  nonequilibrium  population  (note  the  crossing  of  equilibrium  and 
nonequilibrium  curves). 
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3.0  MODELING  RAREFIED  GAS  FLOWS 
3.1  Direct  Boltzmann  Equation  Solution  Methodology 

For  a  continuum  description  of  a  gas  flow,  we  specify  the  velocity  field,  stress  tensor,  and  heat  flux  vector 
and  additionally,  the  three  thermodynamic  variables  of  density,  pressure,  and  temperature.  The  stress 
tensor  and  heat  flux  vector  are  not  usually  treated  as  independent  variables,  such  that  a  minimal  fluid 
mechanical  description  requires  five  independent  variables.  The  statistical  mechanical  description  of  the 
same  system  requires  6N  variables,  where  N  is  the  number  of  molecules  in  the  system;  typical  values  for  N 
are  about  1023.  These  variables  are  usually  taken  to  be  3N  spatial  coordinates,  and  3N  conjugate  momenta 
of  the  constituent  molecules.  As  one  can  see,  the  large  number  of  degrees  of  freedom  in  the  microscopic 
description  distinguishes  it  from  the  macroscopic  description. 

The  Boltzmann  Equation  is  the  governing  equation  in  the  microscopic  description  of  a  dilute  gas  flow,  and 
is  a  nonlinear  integro-differential  equation  for  a  probability  distribution  function  which  describes 
statistically  the  state  of  the  particles  as  a  function  of  time.  The  equation  is  limited  in  its  applicability  to 
systems  for  which  the  following  conditions  are  met: 

•  The  density  is  sufficiently  low  so  that  only  binary  collisions  between  the  constituent  molecules 
need  to  be  considered. 

•  The  spatial  dependence  of  gas  properties  is  sufficiently  low  so  that  collisions  can  be  thought  of  as 
being  localized  in  the  physical  space. 

•  The  intermolecular  potential  is  of  sufficiently  short  range. 

The  Boltzmann  equation  may  be  written  as 


^+PdL 

dt  dx 


'(/) 


(16) 


where 


/(/)  =  LL  (f  f{x4,t)f{x&,t))Bd*d&  (17) 

The  right  side  of  Eq.  (17)  is  a  nonlinear  phase  space  integral  which  describes  the  effect  of  intermolecular 
collisions;  integration  is  performed  over  all  velocity  space  and  over  the  surface  of  a  sphere  of  influence. 
The  collision  term  1(f)  is  a  function  of  the  differential  collision  cross-section  (related  to  B ),  pre-collision 
velocities  (£  £0  and  post-collision  velocities  (£',  £i0-  In  this  representation,  we  assume  both  molecular 
chaos  and  binary  collisions  [24].  A  detailed  explanation  of  the  Boltzmann  collision  term  1(f)  is  beyond  the 
scope  of  the  current  discussion,  but  thorough  descriptions  of  this  term  are  provided  by  Vincenti  and 
Kruger  [8],  Cercignani  [24]  and  Bird  [4].  The  term/,  known  as  the  distribution  function,  describes  the 
many -body  system  and  can  be  integrated  in  velocity  space  to  determine  macroscopic  flow  quantities.  More 
specifically,  f(x,^,t)  is  the  fraction  of  all  particles  in  the  spatial  volume  element  around  x  which  have  their 
velocity  vectors  between  £and  £753£  In  other  words, /is  the  conditional  probability  of  finding  a  particle 
in  the  velocity  space  volume  element  d3^  around  £  provided  that  the  particle  is  located  in  the 
configuration  space  volume  element  d3x  around  x.  Integrating  /  over  its  velocity  argument  will  then  give 
an  expression  for  the  expected  gas  number  density  n  in  the  volume  element  d3x  about  x. 

n{x,t)  =  \x3f(x,^,t)d3^  (18) 


More  generally,  any  macroscopic  flow  quantity  Q  averaged  over  velocity  space  may  be  calculated  as 
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(19) 


where  Q  is  a  function  of  £  and  integrals  in  both  Eqs.  (18)  and  (19)  are  performed  over  all  velocity  space. 
For  example,  the  gas  bulk  velocity  is  found  by  setting  Q  -  £ 

In  a  direct  numerical  solution  of  the  Boltzmann  equation,  Eq.  (16)  is  solved  on  a  grid  with  cell  size  Ax  in 
physical  space  and  on  a  velocity  grid  with  cell  size  introduced  in  a  finite  volume  domain.  On  the  basis 
of  three-dimensional  delta  functions,  the  distribution  function  and  the  collision  integral  can  be  represented 
in  a  discretized  form  [25].  After  determining  the  expansion  coefficients  for  the  collision  integral,  the 
problem  is  reduced  to  solving  the  system  of  equations 


At 


(/; 


(20) 


by  the  finite  difference  method.  Here  y  is  a  phase  space  index,  IY  is  a  discretized  approximation  of  the 
collision  integral,  At  is  the  time  step  interval,  t  is  the  elapsed  simulation  time,  and  A*  indicates  a  difference 
in/between  neighboring  nodes  in  physical  space. 

The  distribution  function  normalized  by  their  respective  peak  values  at  four  x-locations  are  presented  for  a 
Mach  5  shock  wave  in  Figure  6,  as  calculated  using  a  discrete  velocity  Boltzmann  technique.  The  spatial 
locations  are:  far  upstream,  slightly  upstream  of  the  center  of  the  shock,  slightly  downstream  of  the  center 
location,  and  far  downstream.  The  asymptotic  form  of  the  distribution  function  at  the  upstream  and 
downstream  location  is  given  by  the  shifted  Maxwellian  prescribed  by  the  Rankine-Hugoniot  conditions. 
One  can  see  a  non-Maxwellian  distribution  in  x/A  =  0  and  bimodal  distribution  at  x/A  =  0.8.  The 
probability  distribution  function  of  coarse  and  fine  velocity  grids  is  shown  in  Figure  7.  The  relative 
magnitudes  of  the  peaks  are  significantly  different  in  the  shock  wave  (Figure  7a  and  b)  and  close  to  each 
other  far  downstream  (Figure  7c). 

For  hypersonics  applications,  direct  numerical  solution  methods  for  the  discretized  Boltzmann  equation 
have  historically  received  less  attention  than  other  rarefied  flow  simulation  methods,  such  as  direct 
simulation  Monte  Carlo  (DSMC)  [4].  This  can  be  attributed  in  part  to  the  often  prohibitive  computational 
expense  and  memory  requirements  in  solving  the  discretized  Boltzmann  equation  on  a  sufficiently  refined 
velocity  grid.  In  a  recent  study  involving  two  dimensional  simulation  of  a  Mach  4  monatomic  gas  flow 
over  a  cylinder,  a  direct  Boltzmann  calculation  required  roughly  500  times  more  CPU  time  than  a  DSMC 
simulation  [26].  Direct  Boltzmann  methods  do,  however,  have  attributes  which  may  make  them  preferred 
techniques  under  certain  conditions.  In  particular,  these  methods  do  not  suffer  from  the  statistical  scatter 
problems  inherent  in  DSMC,  which  can  make  these  methods  better  suited  to  simulation  of  rarefied 
subsonic  flows  and  to  integration  with  continuum  flow  solvers  in  a  strongly  coupled  hybrid  scheme.  One 
such  hybrid  scheme  is  implemented  in  the  Unified  Flow  Solver  (UFS)  code  [27];  this  code  has  recently 
been  applied  to  a  number  of  hypersonics  problems  involving  mixed  rarefied  and  continuum  flow  regimes, 
and  has  shown  considerable  promise  as  an  engineering  analysis  tool  [28,  29]. 

Several  current  research  topics  for  direct  Boltzmann  methods  are  focused  on  efficiency  improvements. 
These  topics  include  development  of  Boltzmann  schemes  with  higher  order  accuracy  in  velocity  space; 
velocity  grid  adaptation  for  local  optimization  of  the  velocity  grid  based  on  variation  in  Mach  number  and 
temperature;  the  use  of  low  discrepancy  sequences  for  efficient  quasi-Monte  Carlo  evaluation  of  the 
Boltzmann  collision  integral;  and  preferential  sampling  of  high  probability  velocity  space  nodes  during 
collision  integral  evaluation. 
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3.2  Gas  Kinetic  Navier-Stokes  Scheme 

The  Boltzmann  equation  expresses  the  behavior  of  many -particle  kinetic  system  in  terms  of  the  evolution 
equation  for  the  single  particle  gas  distribution  function.  The  simplified  collision  model  given  by  the 
Bhatnagar  Gross  Krook  (BGK)  model  [30]  is  formulated  as 


where /(o)  is  the  equilibrium  (i.e.  Maxwellian)  distribution  function  for  which  integration  over  all  velocity 
space,  as  in  Eq.  (18),  gives  the  local  number  density  n(x,t).  The  symbol  r  in  Eq.  (21)  represents  a 
characteristic  translational  relaxation  time,  and  is  generally  set  to  the  reciprocal  of  the  mean  collision 
frequency  at  position  x  and  time  t.  As  described  above,  the  distribution  function /gives  the  number  density 
per  unit  velocity  space  for  molecules  at  position  x,  velocity  £  and  time  t.  The  left  hand  side  of  Eq.  (21) 
represents  the  free  streaming  of  molecules  in  space,  and  the  right  side  denotes  the  collision  term.  If  the 
distribution  function  /  is  known,  then  macroscopic  variables  of  the  mass,  momentum,  energy  and  stress 
can  be  obtained  by  integration  over  the  moments  of  molecular  velocity.  In  the  BGK  model,  the  collision 
operator  involves  a  simple  relaxation  to  a  local  equilibrium  distribution  function /(0).  The  BGK  model  was 
proposed  to  describe  the  essential  physics  of  molecular  interactions,  while  neglecting  details  of  the 
nonlinear  collision  integral  in  the  full  Boltzmann  equation.  Based  on  the  above  BGK  model,  the  Navier- 
Stokes  equations  can  be  obtained  with  the  Chapman-Enskog  expansion  (discussed  in  Section  3.3) 
truncated  to  the  first  order, 


f  =  fm+Knfm  =  Ao) 


-T 


(  f  f  ^ 

Ao)  [  g  J(Q) 
v  dt  dx  j 


(22) 


where /(o)  and  Knfo  denote  equilibrium  and  nonequilibrium  components,  respectively,  of  the  distribution 
function/.  With  modifications  to  the  characteristic  collision  time  r,  the  validity  of  the  kinetic  model  is 
extended  beyond  that  of  the  traditional  Navier-Stokes  equations  [31,  32]. 

A  kinetic  Navier-Stokes  method  based  on  moments  of  the  BGK  equation,  with  its  approximate  collision 
terms,  is  very  simple  and  attractive.  Attributes  include  the  mathematical  simplicity  of  the  collision  term, 
and  the  lack  of  specific  assumptions  on  intermolecular  forces  of  the  collision  term.  The  BGK  collision 
approximation  is  not  considered  rigorous  theory,  and  one  main  disadvantage  of  the  BGK  method  is  that 
the  relaxation  time  z(x,t)  must  often  be  determined  empirically  for  lack  of  a  better  way. 

Figure  8  shows,  for  a  Mach  3  shock  wave,  the  profiles  of  density,  temperature,  viscous  stress  (xx- 
component)  and  streamwise  component  of  heat  flux.  The  agreement  of  the  gas  kinetic  scheme  results 
agree  well  with  Ref.  [33].  Figure  9  shows,  for  the  Mach  5  shock  wave,  the  density,  mean  temperature,  and 
the  streamwise  and  normal  temperatures  for  the  kinetic  Navier-Stokes  scheme  and  the  Boltzmann  solution. 
From  a  comparison  of  the  density  profiles  in  Figure  9a,  one  notes  the  excellent  agreement  of  the  density 
profiles  for  the  Mach  5  shock  wave.  Note  from  Figure  9b  and  c  that  the  starting  location  of  elevated 
temperatures  predicted  by  the  kinetic  Navier-Stokes  scheme  is  more  upstream  than  that  of  the  Boltzmann 
solution. 


3.3  Moment  Methods 

Another  simplification  to  the  Boltzmann  equation  is  the  well  known  Moment  method  described  briefly  in 
this  paper  from  a  historical  perspective.  Hilbert  [34]  discussed  the  existence  and  uniqueness  of  the 
Boltzmann  equation,  Eq.  (16),  by  exploring  a  restricted  class  of  solutions  (which  can  be  represented  by  a 
power  series,  involving  a  perturbation  parameter),  and  showed  that  there  exists  a  one-to-one  mapping 
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between  the  distribution  function  /  which  is  a  solution  of  the  Boltzmann  equation  and  the  first  five 
moments  of  /  as  given  by  Eq.  (19).  These  moments  represent  density,  bulk  velocity  components  and 
temperature.  Hilbert’s  approach,  although  useful  for  obtaining  Euler  equations  using  the  lowest  order 
approximation,  is  impractical  for  obtaining  the  Navier-Stokes  equations.  In  order  to  derive  both  Euler  and 
Navier-Stokes  equations  from  the  Boltzmann  equation,  an  approach  developed  independently  by  Chapman 
[35,  36]  and  Enskog  [37]  can  be  used.  In  the  Chapman-Enskog  approach,  the  distribution  function /is 
expanded  as  a  power  series  about  the  equilibrium  distribution  function /(0)  as 

f  =  /(o)  +£  f(l)  +£ 2f(2)  + •  •  •  (23) 


where  s  can  be  directly  related  to  the  Knudsen  number.  In  this  approach,  it  is  also  assumed  that  time  is  not 
an  explicit  argument  of  the  distribution  function  /  and  the  time  derivative  of  the  first  five  macroscopic 
moments,  but  is  an  implicit  argument  by  virtue  of  the  dependence  of  /  on  the  macroscopic  moments  and 
their  spatial  derivatives.  This  key  assumption  results  in  expansions  that  are  quite  different  from  those 
obtained  from  Hilbert’s  approach,  especially  in  the  splitting  of  the  time  derivative  off  among  various 
terms  of  the  expansion. 

A  set  of  moment  equations  may  be  derived  by  substituting  the  right  side  of  Eq.  (23)  for /  multiplying  each 
term  of  Eq.  (16)  by  powers  of  <£  then  integrating  over  all  velocity  space.  It  may  be  noted  that  in  the  near 
equilibrium  limit  (as  s  tends  to  zero),  the  Maxwellian  equilibrium  distribution  function  /(o)  is  recovered, 
for  which  moment  equations  result  in  the  Euler  equations,  wherein  the  viscous  stress  tensor  and  the  heat 
flux  vector  are  zero.  Using  a  first  order  approximation  of  the  Chapman-Enskog  expansion  of  the 
Boltzmann  equation,  the  Navier-Stokes  equations  can  be  recovered  from  the  moment  equations,  where 
non-zero  values  of  viscous  stress  tensor  and  heat  flux  vector  (with  dependence  on  the  mean  strain  rate  and 
temperature  gradient  respectively)  are  present.  Second-order  Chapman-Enskog  expansions  applied  to  the 
Boltzmann  equation  will  result  in  moment  equations  called  the  Burnett  equations,  while  higher  order 
expansions  result  in  the  super-Burnett  equations.  The  use  of  Burnett  and  super-Burnett  equations  can  be 
challenging,  as  these  equations  are  ill-posed  and  subject  to  numerical  instabilities. 

One  of  the  models  often  used  to  approximate  the  Boltzmann  collision  term,  Eq.  (17),  is  the  Bhatnagar- 
Gross-Krook  model  [30]  mentioned  in  Section  3.2,  which  greatly  simplifies  the  solution  of  the  Boltzmann 
equation  and  the  derivation  of  the  evolution  of  macroscopic  properties.  Another  interesting  moment  based 
approach  was  proposed  by  Grad  [38],  where  the  distribution  function /is  expanded  in  terms  of  a  complete 
set  of  orthogonal  polynomials  (e.g.  Hermite  polynomials),  whose  coefficients  correspond  to  velocity 
moments.  This  expansion  is  truncated  by  retaining  terms  up  to  a  certain  selected  order,  resulting  in  closed 
set  of  moment  equations  up  to  this  order. 

Although  Grad’s  moment  approach  has  been  widely  used  for  solution  of  the  Boltzmann  equation,  it  is  best 
suited  to  problems  for  which  the  velocity  distribution  function  can  be  expressed  as  a  perturbation  about  the 
equilibrium  (Maxwellian)  distribution.  In  other  words,  Grad’s  method  is  not  very  efficient  for  treating 
problems  with  general  non-equilibrium  (non-Maxwellian)  distributions,  as  Grad’s  moment  method 
assumes  that  the  zeroth  order  term  of  the  expansion  recovers  the  equilibrium  distribution  by  virtue  of  the 
properties  of  Hermite  polynomial  (basis)  functions. 

In  order  to  address  such  problems  for  treatment  of  non-equilibrium  flows,  it  is  plausible  that  an  adaptive 
quadrature  approach  -  where  the  form  of  the  distribution  function  is  not  restricted  to  perturbations  about 
an  equilibrium  distribution  function  -  could  offer  certain  advantages  over  conventional  moment  based 
approaches  (e.g.  Grad’s  moment  approach).  The  implementation  of  such  an  approach  also  requires  a 
detailed  and  efficient  treatment  of  the  collision  operator.  One  such  approach  is  the  application  of  the  direct 
quadrature  method  of  moments  (DQMOM)  for  solution  of  the  full  Boltzmann  equation.  In  DQMOM,  the 
velocity  distribution  function  is  approximated  as  a  collection  of  discrete  Dirac  delta  functions,  with 
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associated  weights  and  abscissas  in  phase  space.  In  particular,  the  particle  distribution  function  is  given  by 
a  discrete  sum  of  N  Dirac  delta  functions,  as 


f(x,^,t)  =  YJPi(x,t)S(4~  £  (x,  t ))  (24) 

(=1 

where  <5(£-<£)  =  S(u  -  Ui)S(v  -  Vj)S(w  -  wj),  ( u ,  v,  w)  are  the  components  of  £  and  the  unknowns  pt  and 
denote  the  associated  quadrature  weights  and  abscissas  in  velocity  space,  respectively.  Each  quadrature 
node  i=  1.  ..  N  has  four  unknowns:  (pi,  uh  vL,  and  hence  there  are  AN  unknowns  to  be  determined  for 
each  spatial  location  x  at  any  point  in  time.  The  unknown  weights  and  abscissas  are  obtained  as  solutions 
of  their  evolution  equations  (or  moment  equations  in  case  of  QMOM),  which  are  in  turn  derived  from  the 
Boltzmann  equation  using  constraints  on  generalized  moments  of  velocity.  These  constraints  on  high  order 
moments  are  derived  based  on  the  complex  collision  integral  given  as  Eq.  (17),  without  resorting  to  the 
use  of  simplifying  kinetic  model  approximations. 

An  example  of  the  generalized  moment  of  velocity  £  with  Cartesian  components  ( u ,  v,  w),  is  given  by 
(upvqw)  where  (p ,  q ,  r)  denote  selected  non-negative  integers.  The  contributions  to  the  evolution  of 
generalized  moments  due  to  collisions  involve  eight  dimensional  integrals.  It  can  be  shown  that  these 
integrals  may  be  evaluated  using  the  above  equation  via  algorithmic  enumeration  of  terms  resulting  from 
products  of  multinomials,  along  with  analytical  evaluation  of  elementary  integrals.  Owing  to  the  adaptive 
and  non-perturbative  nature  of  this  approach  -  unlike  other  approaches  based  on  Knudsen  number 
expansions  -  it  appears  reasonable  to  speculate  that  this  approach  could  be  capable  of  predictions  over  a 
broader  range  of  Knudsen  numbers  with  consideration  of  fewer  moment  constraints. 

The  potential  utility  of  the  adaptive  quadrature  approach  can  be  demonstrated  by  investigation  of 
canonical  problems,  including  homogeneous  relaxation  and  shock  tube  problems.  For  the  case  of  a 
homogeneous  relaxation,  the  behavior  of  the  time  evolution  of  velocity  variances  and  covariances,  shown 
in  Figure  10,  obtained  from  the  quadrature -based  approach  (based  on  the  full  Boltzmann  collision  integral) 
is  found  to  be  in  good  agreement  with  that  obtained  from  DSMC.  Although  not  shown  in  the  figure, 
evolution  of  higher  order  moments  are  also  correctly  captured  using  the  quadrature -based  approach. 
Application  of  the  quadrature  based  approach  to  a  shock  tube  problem,  as  shown  in  Figure  11,  indicates 
that  there  is  a  very  good  agreement  between  the  results  based  on  the  Boltzmann  collision  operator  and  its 
BGK  approximation  at  low  Knudsen  numbers. 

3.4  Direct  Simulation  Monte  Carlo 

The  direct  simulation  Monte  Carlo  (DSMC)  method  is  the  most  widely  used  technique  for  numerical 
modeling  of  rarefied  gas  flows,  and  over  the  past  several  decades  DSMC  has  been  an  important  tool  in 
design  and  analysis  for  a  wide  variety  of  hypersonic  flow  problems  [4].  The  DSMC  method  was  initially 
developed  by  Bird  in  the  early  1960’s  for  use  in  simple  homogeneous  relaxation  problems,  in  which  a 
nonequilibrium  gas  velocity  distribution  is  driven  toward  equilibrium  by  intermolecular  collisions  [39].  In 
the  decades  since,  DSMC  has  been  applied  with  great  success  to  a  wide  variety  of  hypersonic  engineering 
problems  [40,  41,  42]. 

While  the  underlying  assumptions  in  DSMC  are  valid  for  any  dilute  gas  flow  (i.e.  a  gas  flow  for  which  the 
collision  cross  section  is  much  smaller  than  the  spacing  between  atoms  or  molecules)  the  application  of 
DSMC  has  generally  been  limited  to  rarefied  flows  due  to  the  higher  efficiency  of  continuum  CFD 
techniques  for  continuum  flow  simulation.  Within  the  transition  Knudsen  number  regime,  DSMC 
combines  a  high  degree  of  accuracy,  conceptual  simplicity,  ease  of  application  to  complex  geometries,  and 
unconditional  numerical  stability.  For  high  speed  rarefied  flows,  DSMC  also  tends  to  be  considerably  less 
expensive  than  other  rarefied  flow  simulation  methods,  and  the  applicability  of  DSMC  to  complex  three  - 


RTO-EN-AVT-194 


1  -13 


Review  of  Rarefied  Gas  Effects  in  Hypersonic  Applications 


ORGANIZATION 


dimensional  flows  has  increased  over  time  due  to  continual  improvements  in  available  computer 
resources.  Along  with  development  and  availability  of  full-featured  DSMC  codes  [43-47],  and  particularly 
the  freely  available  codes  of  Bird  [48],  these  attributes  tend  to  make  DSMC  a  preferred  tool  for  use  in  both 
design  and  aerothermal  analysis  of  hypersonic  vehicles. 

3.4.1  Overview  of  Numerical  Procedures 

At  its  core,  the  DSMC  method  involves  tracking  and  collision  operations  among  a  collection  of  simulated 
particles,  each  representing  a  large  number  of  atoms  or  molecules,  in  a  manner  which  reproduces  the 
underlying  physics  of  the  governing  Boltzmann  equation  [4].  In  a  typical  DSMC  simulation,  every  particle 
is  assigned  a  position,  a  velocity,  a  number  or  pointer  to  identify  the  particle  species,  and  internal 
(rotational  and  vibrational)  energy  values.  Particles  are  grouped  into  finite  volume  computational  cells  for 
use  in  both  collision  partner  selection  and  time-averaged  sampling  of  simulation  results.  Within  each  cell, 
the  gas  density  is  calculated  as  a  function  of  the  number  of  particles  per  unit  volume,  the  gas  bulk  velocity 
is  the  mean  velocity  of  all  particles  in  that  cell,  and  the  translational  temperature  is  proportional  to  the 
mean-square  particle  speed  in  a  coordinate  frame  which  moves  at  the  bulk  velocity.  Rotational  and 
vibrational  temperatures  -  which  may  differ  significantly  from  the  translational  temperature  in  a  rarefied 
hypersonic  flow  -  are  computed  from  the  mean  internal  energy  of  all  particles  in  the  cell.  Surface  fluxes 
are  calculated  from  the  momentum  and  energy  exchange  which  occurs  as  particles  are  reflected  off  of 
solid  wall  boundaries. 

DSMC  is  an  inherently  explicit  technique,  and  a  DSMC  simulation  of  a  steady  state  flow  typically  begins 
with  either  vacuum  or  uniform  initial  conditions.  The  flowfield  then  evolves  over  time  toward  steady  state, 
after  which  time-averaged  sampling  may  be  performed  over  a  number  of  time  steps  to  reach  desired  levels 
of  statistical  scatter  in  the  output  quantities  of  interest. 

During  each  simulation  time  step,  a  typical  DSMC  algorithm  involves  three  basic  procedures:  First, 
particles  are  moved  through  a  computational  grid  according  to  assigned  velocities,  and  are  reassigned  as 
necessary  to  new  cells  in  this  grid.  Particles  which  exit  the  simulation  domain  through  inflow  or  outflow 
boundaries  are  removed,  while  any  particle  which  collides  with  a  solid  wall  boundary  is  reflected  off  this 
boundary  according  to  some  gas-surface  interaction  model.  Additional  particles  may  then  be  introduced 
along  inflow  boundaries.  Second,  particles  within  each  cell  are  organized  into  pairs,  and  collision 
operations  are  performed  on  these  pairs  in  a  manner  which  produces  both  the  desired  collision  frequency 
and  desired  macroscopic  transport  coefficients.  Energy  exchange  between  translational  and  internal  energy 
modes,  as  well  as  chemical  reactions,  may  be  applied  to  some  fraction  of  colliding  particles.  Finally,  if 
either  steady  state  conditions  have  been  reached  or  an  unsteady  flow  is  being  simulated,  then  sampling 
routines  are  used  to  compile  relevant  field  and  surface  quantities  for  output  of  simulation  results. 

These  basic  DSMC  procedures,  which  are  performed  for  every  cell  during  each  time  step,  are  outlined  in 
Figure  12.  Note  that  the  order  of  these  procedures  varies  between  different  implementations  of  DSMC, 
and  the  detailed  algorithms  and  physical  models  associated  with  each  procedure  can  vary  considerably  as 
well. 


3.4.2  Conditions  for  Accurate  Simulation 

For  accurate  DSMC  simulation,  two  numerical  criteria  must  usually  be  met.  First,  the  time  step  interval 
must  be  small  compared  to  the  mean  collision  time.  Second,  the  mean  distance  between  colliding  particles 
-  often  termed  the  “mean  collision  separation”  (MCS)  -  should  be  small  in  relation  to  the  local  mean  free 
path  (MFP)  [4,  48].  The  first  criterion  follows  from  an  operator-splitting  approach  used  to  decouple 
particle  movement  and  collisions  during  each  time  step;  this  criterion  is  met  either  by  careful  limiting  of  a 
global  time  step  interval,  or  by  local  time  step  adaptation  with  different  time  step  intervals  employed  in 
different  cells.  The  second  criterion,  based  on  a  desired  MCS/MFP  ratio,  is  often  more  difficult  to  meet 
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without  compromising  simulation  efficiency,  and  has  been  the  subject  of  considerable  research  over  the 
past  few  decades.  Two  basic  strategies  are  generally  utilized  (often  together)  to  limit  the  MCS  in  a  DSMC 
simulation:  grid  adaptation,  and  reducing  grid  dependence  in  the  MCS  through  either  adaptive  subcells 
[48]  or  nearest  neighbor  collisions  [49].  Both  strategies  are  discussed  below. 

Variation  in  grid  architecture  and  grid  adaptation  techniques  are  arguably  the  greatest  distinguishing 
characteristics  among  the  different  widely  used  DSMC  codes.  Unstructured  body-fitted  grids,  as  employed 
in  the  MONACO  DSMC  code  [44],  provide  the  highest  degree  of  flexibility,  but  require  potentially  time- 
consuming  grid  generation  by  the  user  and  can  complicate  procedures  for  automatic  adaptation  to  the  local 
mean  free  path.  Most  other  codes  utilize  some  form  of  Cartesian  grid,  with  cuboidal  or  rectangular  cells, 
cut  cells  along  inclined  surface  boundaries,  and  capabilities  to  automatically  divide  each  cell  into  a 
number  of  smaller  cells  using  refinement  criteria  based  on  the  local  mean  free  path  [43,  45,  46].  This 
strategy  avoids  requirements  to  import  an  externally  defined  grid,  and  tends  to  ease  simulation  setup  for 
flow  problems  involving  complex  surface  geometries.  Cartesian  grid  architectures  can  be  further  divided 
into  two-level  or  three-level  Cartesian  grids,  tree  structures  [50],  and  multi-element  groupings.  Other  grid 
types  [4],  such  as  a  monotonic  Lagrangian  grid  [51],  have  also  been  used  with  some  success.  Detailed 
descriptions  of  the  various  grid  adaptation  strategies  are  provided  in  the  references. 

While  grid  adaptation  based  on  the  local  mean  free  path  is  the  most  common  way  to  ensure  sufficiently 
small  MCS,  the  MCS  can  be  further  reduced  by  means  of  subcells  or  nearest  neighbor  selection  routines. 
In  transient  adaptive  subcell  procedures,  all  particles  in  a  given  cell  are  organized  into  smaller  subcells  as 
a  first  step  in  collision  procedures  during  each  time  step  [48].  Subcell  sizes  are  based  on  either  the  local 
mean  free  path  or  the  desired  number  of  particles  per  subcell.  Collision  partners  are  then  selected  from  the 
same  subcell  if  possible;  otherwise  particles  in  nearby  subcells  are  paired  together.  In  nearest  neighbor 
selection,  the  first  particle  in  each  collision  pair  is  randomly  selected,  and  the  closest  particle  among  all 
particles  in  the  cell  is  chosen  as  the  second  particle  in  the  pair  [49].  Typical  modifications  to  this 
procedure  involve  preventing  the  same  two  particles  from  colliding  repeatedly,  and  limiting  the  nearest 
neighbor  search  to  some  randomly  selected  subset  of  all  particles  in  the  cell  [53,  54]. 

3.4.3  Collision  Procedures 

As  outlined  in  Figure  12,  DSMC  collision  procedures  include  several  different  steps,  and  an  overview  of 
these  steps  is  provided  in  the  following  discussion.  In  the  first  step  for  collision  procedures  applied  to  a 
given  cell,  some  number  Npairs  of  particle  pairs  are  selected  from  among  all  particles  in  the  cell.  The  value 
of  Npairs  varies  among  different  collision  selection  schemes.  In  the  popular  “no  time  counter”  (NTC) 
scheme  of  Bird  [4],  Npairs  is  calculated  as 


Nmirs=-Nn(cjc)  At  (25) 

where  N  is  the  number  of  particles  in  the  cell,  n  is  the  time-averaged  number  density,  (crc)max  is  the 
maximum  product  of  collision  cross  section  and  relative  speed  for  any  collision  in  the  cell  over  a  large 
number  of  time  steps,  and  At  is  the  time  step  interval.  Each  pair  is  then  considered  to  participate  in  a 
collision;  using  the  NTC  scheme,  a  pair  is  selected  if  cfyc;/(crc)max  >  R ,  where  o;,  and  cLJ  are  respectively  the 
collision  cross  section  and  relative  speed  for  the  (/,  j)  particle  pair,  and  R  is  a  random  number  in  [0,1]. 
Other  schemes,  such  as  the  majorant  frequency  scheme  (MFS)  [45]  and  a  scheme  of  Baganoff  and 
McDonald  [52],  involve  alternate  procedures  for  selecting  particles  to  collide. 

In  NTC  as  in  other  collision  selection  schemes,  collision  probabilities  are  proportional  to  a  collision  cross 
section  cr,  which  is  in  turn  generally  a  function  of  the  relative  incident  speed  c  for  the  two  particles.  The 
dependence  of  cron  c  is  directly  related  to  the  dependence  of  transport  coefficients  -  such  as  viscosity  -  on 
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temperature.  Most  widely  used  DSMC  codes  employ  the  variable  hard  sphere  model  [4],  for  which  a  is  a 
power  law  function  of  c  and  viscosity  scales  with  a  constant  power  of  translational  temperature. 

Once  colliding  pairs  have  been  determined,  each  pair  is  considered  for  probabilistic  redistribution  of 
energy  between  translational  and  internal  energy  modes.  For  any  collision  involving  diatomic  or 
polyatomic  molecules,  some  energy  may  be  exchanged  between  translational  and  rotational  modes,  and 
the  probability  of  such  exchange  is  equal  to  the  reciprocal  of  the  rotational  collision  number.  This  collision 
number  is  either  set  to  some  constant  or  is  allowed  to  vary  as  a  function  of  collision  energy  or  cell- 
averaged  temperature  [55].  If  translational -rotational  energy  exchange  is  determined  to  have  occured,  then 
new  particle  rotational  energies  are  typically  resampled  from  continuous  equilibrium  distributions  using 
the  technique  of  Borgnakke  and  Larsen  [56].  For  high  temperature  flow  simulations  involving  vibrational 
excitation,  similar  probabilistic  procedures  are  employed  to  determine  whether  translational -vibrational 
energy  exchange  occurs,  and  to  resample  vibrational  energy  values  from  quantized  equilibrium 
distributions  [4].  For  application  to  very  high  temperature  flows  or  flows  involving  gas  radiative  emission, 
a  similar  treatment  has  been  applied  to  nonequilibrium  electronic  excitation  [4,  57]. 

Note  that  commonly  used  energy  exchange  procedures  [56]  are  phenomenological,  and  do  not  account  for 
potentially  important  effects  such  as  vibrational-vibrational  exchange,  anharmonic  vibrational  quantum 
levels,  and  state-dependent  exchange  probabilities.  Although  not  usually  included  in  DSMC  hypersonic 
flow  simulations,  all  of  these  phenomena  have  been  considered  in  a  research  context  [58].  It  should  be 
emphasized  that  the  particle  nature  of  DSMC,  with  individual  internal  energy  values  assigned  to  each 
particle,  provides  considerable  potential  for  accurate  and  efficient  simulation  of  gas  flows  involving 
nonequilibrium  rotational,  vibrational  or  electronic  energy  distributions. 

For  accurate  simulation  of  high  enthalpy  hypersonic  flows,  chemistry  modeling  capabilities  are  often 
required.  Several  DSMC  chemistry  models  have  been  developed,  and  of  these  the  total  collision  energy 
and  vibrationally  favored  dissociation  models  are  likely  in  most  use  [4,  58].  These  models  are  designed  to 
reproduce  the  temperature  dependence  of  modified  Arrhenius  reaction  rate  coefficients,  by  calculating  the 
reaction  probability  for  a  given  collision  pair  as  a  function  of  relative  translational  and  internal  energies  of 
the  colliding  particles.  Due  to  the  low  density  in  most  DSMC  applications,  chemistry  modeling  efforts 
have  primarily  focused  on  dissociation  and  exchange  reactions,  although  some  DSMC  codes  include 
capabilities  for  recombination  and/or  surface  catalysis  [4].  Recent  work  of  Bird  and  others  has  focused  on 
integrated  treatment  of  reactions  and  vibrational  state  transitions,  with  the  potential  to  avoid  reliance  on 
near-equilibrium  assumptions  such  as  those  used  in  the  modified  Arrhenius  rates  [59]. 

Following  any  translational-internal  energy  exchange  and  chemistry  procedures  for  a  given  collision  pair, 
the  relative  translational  energy  of  the  two  particles  is  adjusted  in  order  to  enforce  conservation  for  the 
total  energy  in  the  collision.  Finally,  the  translational  energy  is  distributed  among  the  two  particles  in  a 
manner  which  enforces  both  energy  and  momentum  conservation  [4].  For  the  standard  variable  hard 
sphere  (VHS)  collision  models  in  a  simple  (one  species)  gas,  this  involves  randomly  sampling  post¬ 
collision  particle  velocities  from  opposite  points  on  a  spherical  shell  in  velocity  space,  which  is  centered  at 
the  average  incident  velocity  of  the  two  particles.  Other  collision  models  employ  anisotropic  resampling  of 
the  post-collision  relative  velocity,  in  order  to  more  accurately  simulate  transport  phenomena  or  other 
effects  [4].  In  particular,  the  variable  soft  sphere  (VSS)  model  [60]  corrects  for  the  inability  of  the  VHS 
model  to  independently  achieve  desired  viscosity  and  mass  diffusion  rates,  and  the  variable  sphere  (VS) 
[61]  model  simulates  the  effects  of  various  intermolecular  potential  functions  by  fixing  the  deflection 
angle  between  incident  and  post-collision  relative  velocities. 

3.4.4  Current  Research  Issues 

The  DSMC  method  has  advanced  considerably  in  the  past  several  decades,  and  numerous  researchers  have 
made  important  improvements  in  accuracy,  efficiency  and  applicability  of  the  method.  However,  several 
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research  issues  remain.  Due  to  advances  in  computing  power,  DSMC  continues  to  increase  in  popularity 
and  become  more  of  a  practical  tool  for  complex  engineering  problems.  With  the  increasing  value  of 
DSMC  as  an  engineering  tool  for  hypersonic  vehicle  design  and  analysis,  a  number  of  research  problems 
have  recently  received  considerable  attention,  and  current  DSMC  hypersonics  research  is  focused  in 
several  different  areas.  Some  of  these  research  areas  are  cited  in  the  following  discussion. 

One  important  set  of  research  topics  involves  improved  characterization  of  nonequilibrium  internal  energy 
distributions,  with  a  particular  emphasis  on  quantized  vibrational  excitation.  In  a  typical  high-enthalpy 
hypersonic  flow,  highly  nonequilibrium  vibrational  energy  distributions  are  expected  throughout  the  shock 
layer  and  wake,  and  surface  heating,  chemical  composition  and  radiative  properties  may  be  strong 
functions  of  local  vibrational  energy  distributions  within  these  regions.  While  the  traditional  DSMC 
approach  [56,  48]  to  vibrational  excitation  allows  for  arbitrary  energy  distributions,  this  approach  is 
phenomenological  and  neglects  much  of  the  underlying  physical  phenomena.  Standard  DSMC  procedures 
for  vibrational  excitation  involve  probabilistic  sampling  from  equilibrium  distributions  [56],  employ  a 
simple  harmonic  oscillator  approximation,  and  do  not  consider  potentially  important  processes  such  as 
vibrational-vibrational  exchange.  Vibration-dissociation  coupling  is  also  either  neglected  [4]  or  greatly 
simplified  [58]. 

To  correct  for  limitations  in  traditional  DSMC  approaches  for  internal  energy  exchange,  much  interest  has 
recently  focused  on  development  of  state-to-state  quantum  kinetics  models  [62],  and  in  new  chemistry 
models  which  account  for  the  strongly  coupled  nature  of  vibrational  excitation  and  chemical  reaction 
mechanisms  [59].  Figure  13  illustrates  one  example,  taken  from  Ref.  [62],  of  the  potential  accuracy 
improvements  associated  with  a  state-to-state  model  for  vibrational  quantum  level  transitions.  The  figure 
shows  relative  populations  at  quantum  level  4,  as  a  function  of  transit  time  through  a  Mach  7  shock  in  CO. 
Results  from  a  traditional  DSMC  simulation  are  on  the  left,  while  results  from  a  DSMC  simulation 
utilizing  a  new  state-to-state  model  on  the  right.  (Both  DSMC  results  are  shown  in  red.)  In  comparing  with 
experimental  data  points  also  displayed  in  Figure  13,  we  find  a  substantial  improvement  in  accuracy  for 
the  new  model  during  the  nonequilibrium  relaxation  period  at  small  elapsed  time  values. 

Flowfield-radiation  coupling  is  another  active  area  of  DSMC  research,  and  effects  of  shock  layer  radiation 
on  high  enthalpy  reentry  flowfields  may  need  to  be  considered  for  accurate  gas  flow  simulation. 
Potentially  important  coupling  effects  in  such  flows  include  spectrally  resolved  gas  emission  and 
absorption,  as  well  as  radiation  contributions  from  the  vehicle  surface  and  solid  particles  associated  with 
surface  ablation  [63,  64].  Radiative  emission  has  been  linked  to  vibrational  energy  level  transitions,  and 
improved  characterization  of  these  transitions  should  allow  improved  capabilities  for  estimating  both 
radiative  heating  for  reentry  vehicles  and  plume  radiation  signatures. 

Multiphase  rarefied  flows  are  of  interest  in  simulating  high  altitude  solid  rocket  plumes  and  spacecraft 
venting  flows,  and  recent  work  has  focused  on  development  of  DSMC  models  for  two-way  coupled 
interaction  between  a  gas  and  condensed  phase  particles  [65].  Related  work  has  involved  model 
development  for  nucleation,  condensation  and  evaporation  of  molecular  clusters  in  freely  expanding 
plumes. 

As  mentioned  above,  several  different  strategies  have  been  employed  for  DSMC  grid  adaptation  [4],  and 
recent  work  on  multi-level  Cartesian  grid  refinement  has  shown  great  promise  in  effectively  adapting  grids 
for  complex  flowfields  [46].  Another  newly  proposed  strategy  employs  a  uniform  Cartesian  grid  with 
modified  collision  procedures  to  reduce  cell  size  dependence  in  simulation  results  [66].  The  DSMC 
community  is  far  from  a  consensus  on  the  most  effective  grid  adaptation  technique,  and  this  remains  an 
ongoing  research  topic. 

For  efficient  and  accurate  simulation  of  hypersonic  flows  involving  a  wide  range  of  Knudsen  number 
regimes,  much  recent  work  has  involved  development  of  hybrid  schemes  which  utilize  DSMC  only  in 
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rarefied  regions  and  apply  Navier-Stokes  techniques  [67,  68,  69]  or  alternate  particle  methods  [70,  71] 
over  continuum  portions  of  the  flowfield.  Such  hybrid  approaches  allow  DSMC  to  be  employed  only 
within  regions  where  DSMC  is  relatively  efficient  and  where  approximations  underlying  continuum 
approaches  are  assumed  to  be  invalid.  Hybrid  techniques  are  particularly  useful  for  simulating  low  altitude 
hypersonic  flows,  where  full  DSMC  simulation  may  be  prohibitively  expense  but  strong  local 
nonequilibrium  occurs  in  the  velocity  distribution  within  small  flowfield  regions.  These  regions  include 
the  bow  shock,  portions  of  the  boundary  layer  over  leeward  surfaces,  the  area  surrounding  a  sharp  leading 
edge,  and  low  density  portions  of  the  wake. 

Particular  challenges  in  development  of  DSMC-based  hybrid  schemes  for  simulation  of  this  type  of 
“multi-scale”  flow  include  accurate  determination  of  continuum  breakdown  for  domain  decomposition 
between  DSMC  and  continuum  modules  of  a  hybrid  code,  and  communication  between  the  two  modules 
along  the  continuum  breakdown  boundary.  The  gradient  length  Knudsen  number  Knmax ,  given  by  Eq.  (7), 
is  often  employed  to  automatically  divide  the  computational  domain  between  DSMC  and  continuum 
regions,  with  some  cutoff  value  (typically  0.05)  used  to  indicate  continuum  breakdown.  However,  accurate 
assessment  of  continuum  breakdown  can  require  alternate  breakdown  parameters,  and  the  optimal  cutoff 
value  may  vary  widely  between  different  flows.  Several  alternate  breakdown  parameters  have  been 
proposed,  and  characterization  of  continuum  breakdown  is  currently  an  active  research  topic. 
Communication  between  DSMC  and  continuum  methods  is  another  area  of  current  research,  and  is  greatly 
complicated  by  the  potentially  large  statistical  scatter  in  conserved  macroscopic  quantities  as  calculated  by 
DSMC.  Approaches  to  compensate  for  this  scatter  include  the  use  of  additional  DSMC  particles  within 
“ghost  cells”,  time -averaging  for  DSMC  quantities  [67,  72],  and  continuum  particle  methods  which  are 
less  sensitive  to  scatter  effects  than  typical  Navier-Stokes  schemes  [70,  71]. 

Another  current  DSMC  research  area  of  particular  importance  to  rarefied  hypersonic  flows  is  gas -surface 
interaction.  The  standard  Maxwell  model  for  such  interactions,  which  is  implemented  in  most  DSMC 
codes,  uses  a  constant  thermal  accommodation  coefficient  (TAC)  to  specify  the  probability  that  any 
particle-wall  boundary  collision  results  in  diffuse  reflection;  otherwise  the  collision  will  involve  specular 
reflection  [4].  A  diffusely  reflected  particle  is  assigned  new  velocity  components  (and  typically  new 
internal  energy  values  as  well)  which  are  sampled  from  distributions  corresponding  to  a  stationary  gas 
reservoir  at  the  wall  boundary  temperature.  For  a  specularly  reflected  particle,  the  normal  velocity 
component  is  reversed  while  tangential  velocity  components  are  unchanged.  While  the  Maxwell  model 
has  positive  attributes  such  as  simplicity  of  implementation  and  availability  of  TAC  values  for  several 
surface-gas  combinations,  the  model  is  very  approximate  and  potentially  inaccurate.  In  particular,  this 
model  neglects  complex  scattering  characteristics  observed  through  molecular  beam  experiments  and 
detailed  molecular  dynamics  simulations.  Alternate  models,  such  as  the  Cercignani-Lampis-Lord  (CLL) 
model  [73]  add  physical  realism  through  the  use  of  additional  parameters,  but  have  yet  to  gain  widespread 
use  due  to  a  general  lack  of  knowledge  regarding  appropriate  input  parameter  values.  For  realistic 
treatment  of  gas-surface  interaction,  an  ideal  model  would  meet  several  different  criteria.  These  include 
enforcement  of  the  reciprocity  relation  [4];  consideration  of  surface  roughness  and  temperature  effects; 
depends  on  gas  species  and  surface  composition;  independent  accommodation  of  internal  energy  modes; 
surface  catalysis,  sublimation  and  deposition;  and  scattering  distributions  which  agree  well  with 
experimental  data.  As  assumptions  and  input  parameters  related  to  gas-surface  interaction  may  have 
considerable  effects  on  hypersonic  vehicle  characteristics  -  particularly  surface  heating  and  aerodynamic 
forces  -  additional  work  in  both  DSMC  model  development  and  experiment  seems  warranted. 

Several  other  active  research  areas  exist  within  the  DSMC  community.  These  include  development  of 
DSMC-based  procedures  for  efficient  simulation  of  subsonic  flows  [74];  alternate  operator  splitting 
schemes  for  particle  advection  and  collision  processes  [75];  efficient  minimization  of  the  mean  collision 
spacing  [53,  54];  analysis  of  temporal  and  spatial  discretization  errors  [76];  ionization,  electronic 
excitation  and  electron  transport  using  either  particle  or  finite  volume  techniques  [77,  78,  79];  time  step 
and  numerical  weight  adaptation  [4];  improved  modeling  of  collision  dynamics  at  very  low  or  very  high 
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temperatures  [80];  identification  of  convergence  to  steady  state  [81];  and  techniques  to  better  utilize 
computer  resources,  such  as  DSMC  implementations  for  graphics  processing  units,  combined  use  of 
shared  and  distributed  memory,  and  improved  procedures  for  parallel  domain  decomposition  [82]. 

3.4.5  DSMC  Simulation  Examples 

Over  the  past  several  decades,  the  DSMC  method  has  been  used  extensively  in  a  variety  of  hypersonic 
applications.  An  exhaustive  list  of  these  applications  is  out  of  the  scope  of  this  work,  but  a  number  of 
representative  examples  are  mentioned.  These  examples  are  organized  under  the  general  categories  of 
aerothermal  analysis  for  atmospheric  entry  and  hypersonic  cruise  vehicles;  and  thrusters,  plumes  and 
spacecraft  aerodynamics.  Flows  in  the  first  category  are  distinguished  by  relatively  large  aerodynamic 
drag  and  heating  effects,  and  vehicle  geometries  which  are  primarily  designed  with  aerothermal 
characteristics  in  mind.  Flows  in  the  second  category  tend  to  involve  near-free  molecular  freestream 
conditions,  with  the  highest  gas  densities  occurring  in  nearfield  plume  regions  for  thrusters  or  venting 
sources. 

3.4.5. 1  Aerothermal  Analysis  for  Atmospheric  Entry  and  Hypersonic  Cruise  Vehicles 

Due  to  the  difficulty  and  expense  of  flight  experiments  involving  high  Knudsen  number  flow  regimes, 
most  experimental  comparisons  for  assessment  of  DSMC  codes  and  techniques  are  based  on  subscale 
wind  tunnel  data.  Likely  the  most  comprehensive  of  these  comparisons  was  promoted  by  Working  Group 
18  of  the  NATO  AGARD  Fluid  Dynamics  Panel,  and  involved  wind  tunnel  tests  at  four  different 
experimental  facilities  and  extensive  comparison  among  several  DSMC  and  Navier-Stokes  codes.  A 
spherically  blunted  70°  cone,  with  the  same  forebody  proportions  as  the  Mars  Pathfinder  aeroshell,  was 
tested  under  a  variety  of  conditions  involving  Mach  numbers  between  about  10  and  20.  Most  test  cases 
used  N2  gas,  with  sufficiently  low  freestream  enthalpy  for  negligible  dissociation  and  relatively  small 
vibrational  excitation  effects.  Measured  data  included  surface  pressure  and  heat  flux,  gas  density 
throughout  the  flowfield,  and  normal  and  axial  forces  for  a  range  of  pitch  angles.  Rarefaction  effects  were 
observed  in  both  the  wake  and  the  forebody  shock  layer,  and  good  overall  agreement  was  generally  found 
between  experimental  data  and  DSMC.  Study  results  are  reviewed  in  Ref.  [83]. 

Two  examples  of  comparisons  between  DSMC  and  experimental  data  from  the  NATO-AGARD  study,  as 
taken  from  Ref.  [41],  are  shown  in  Figure  14.  Contours  of  normalized  density  are  displayed  from 
axisymmetric  MONACO  [44]  DSMC  results  and  from  electron  fluorescence  measurements  of  Allegre  et 
al.  [84]  for  a  Mach  20.2  case  at  Kn  =  0.03  and  0°  angle  of  attack.  Centerline  heat  flux  results  are  compared 
for  a  somewhat  higher  enthalpy  condition  at  Mach  15.6,  Kn  =  0.002  and  0°  angle  of  attack.  Figure  15 
shows  pressure  and  density  contours  for  a  case  at  Mach  20.2,  Kn  =  0.03  and  an  angle  of  attack  of  10°.  The 
pressure  plot  includes  surface  pressure  and  scalar  field  pressure  contours  as  well  as  streamlines,  and  was 
generated  from  a  three  dimensional  simulation  using  the  newly  developed  HAP  DSMC  code  [66]. 
Experimental  number  density  values  in  Figure  15  are  normalized  by  the  local  density  in  an  empty  test 
section.  Significant  asymmetry  is  observed  in  both  forebody  and  wake  regions  due  to  the  angle  of  attack, 
and  relatively  good  agreement  is  found  between  DSMC  and  experiment. 

In  another  comprehensive  comparison  of  DSMC  and  ground  test  data,  a  series  of  axisymmetric  hypersonic 
flows  over  double  cone  and  cylinder-flare  configurations  were  used  to  investigate  shock-shock  and  shock¬ 
boundary  layer  interactions  under  rarefied  conditions  [85].  Representative  results  from  this  study  are 
shown  in  Figure  16,  for  a  case  involving  Mach  15.6  flow  of  nitrogen  over  a  25755°  double  cone  with 
Kn  =  0.002.  A  plot  of  pressure  contours  and  streamlines  shows  extremely  complex  shock-shock 
interactions  near  the  junction  between  the  two  cones,  with  a  clearly  defined  recirculating  region.  Also  in 
Figure  16,  experimental  data  for  surface  heat  flux  is  compared  with  simulation  results  from  two  different 
DSMC  codes,  and  very  good  overall  agreement  is  found.  Note  the  reduction  in  heat  flux  around  the 
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recirculation  zone,  and  a  heat  flux  spike  where  a  secondary  shock  intersects  the  boundary  layer. 

As  an  example  of  DSMC  results  for  a  realistic  three  dimensional  flow  over  a  full  scale  hypersonic  vehicle, 
Figure  17  shows  contours  of  translational  temperature,  surface  heat  flux,  Mach  number,  and  streamlines 
for  a  Mach  14  flow  of  air  over  a  viscous-optimized  waverider  geometry  [86]  at  0°  angle  of  attack  and  a 
simulated  altitude  of  100  km.  The  global  Knudsen  number,  based  on  a  body  length  of  5  m,  is  0.02.  This 
simulation  was  performed  using  the  HAP  DSMC  code  [66],  in  order  to  evaluate  degradation  in  waverider 
lift  expected  at  very  high  altitudes  due  to  viscous  and  rarefaction  effects.  The  computed  lift-to-drag  ratio 
of  about  0.3  was  well  below  that  expected  at  lower  altitudes. 

Figure  18  shows  streamlines  and  pressure  contours  from  a  DSMC  simulation  of  an  inflatable  “ballute” 
concept  to  increase  high  altitude  deceleration  during  atmospheric  reentry,  as  taken  from  Moss  [87]. 
Complex  shock-shock  interactions  are  found  to  result  from  shocks  generated  around  opposite  sides  of  the 
toroidal  ballute,  and  around  the  capsule  to  which  the  ballute  is  tethered.  In  a  detailed  DSMC  study  of  high 
altitude  ballute  aerothermodynamics  [87],  a  variety  of  flowfield  and  surface  flux  characteristics  were 
assessed. 

Mach  number  contours  for  a  Mach  10  flow  of  nitrogen  over  a  subscale  Apollo  command  module  geometry 
are  shown  in  Figure  19,  as  taken  from  DSMC  results  used  in  a  comparison  with  aerodynamic  coefficients 
from  a  series  of  wind  tunnel  measurements  [88].  Here  the  Knudsen  number  is  0.067  and  the  angle  of 
attack  is  30°.  Significant  discrepancies  were  observed  between  DSMC  and  experimental  data  for  both 
integrated  vehicle  drag  and  pitching  moment,  although  relatively  good  agreement  was  found  for  a  slender 
cone  geometry.  Excellent  agreement  was  also  found  between  aerodynamic  coefficients  in  Ref.  [88]  and 
those  calculated  using  another  DSMC  code,  as  discussed  in  Ref.  [41]. 

Figure  20  shows  the  surface  pressure  distribution  at  an  altitude  of  400,000  ft  for  the  X-38  vehicle  as 
computed  using  the  DAC  DSMC  code  [40,  43].  The  X-38  was  designed  as  a  crew  return  vehicle  for  the 
International  Space  Station,  and  DSMC  simulations  were  used  to  evaluate  aerodynamic  characteristics  for 
reentry  conditions  between  free  molecular  and  continuum  flow  regimes. 

Contours  of  normalized  density  for  the  Space  Shuttle  Orbiter  at  170  km  altitude  are  shown  in  Fig.  21,  as 
taken  from  DSMC  results  of  Rault  [89].  DSMC  simulations  were  performed  for  a  range  of  shuttle  reentry 
conditions  between  100  km  and  170  km  altitudes,  and  DSMC  aerodynamic  coefficients  were  compared 
with  flight  data,  wind  tunnel  data,  and  with  Navier-Stokes  and  bridging  formula  calculations. 

In  a  recent  example  of  the  utility  of  DSMC  for  hypersonic  stability  analysis,  Abdel -jawad  et  al.  [90] 
performed  a  series  of  DSMC  simulations  for  the  ESA  Beagle2  vehicle  during  Mars  atmospheric  entry,  and 
used  simulation  results  within  a  vehicle  dynamics  model  to  demonstrate  static  instability.  The  authors 
conclude  that  simple  bridging  relations,  which  had  been  employed  in  pre-flight  analysis  to  determine  the 
required  rate  of  spin  during  entry,  had  underestimated  the  effect  of  viscous  forces  and  may  have  ultimately 
led  to  mission  failure.  Figure  22  shows  overall  temperature  contours  from  a  DSMC  simulation  used  in  this 
study,  at  an  angle  of  attack  of  1 1°  and  Kn  =  0.0291. 

3  A.  5. 2  Thrusters,  Plumes  and  Spacecraft  Aerodynamics 

In  addition  to  hypersonics  applications  involving  high  altitude  aircraft  and  reentry  vehicles,  common 
DSMC  applications  include  exo-atmospheric  plume  impingement  analysis  and  satellite  aerodynamics.  One 
such  example  is  the  extensive  use  of  the  SMILE  DSMC  code  [45]  for  transition  regime  analysis  of  the  Mir 
space  station  during  atmospheric  descent,  prior  to  its  controlled  deorbit  in  2001  [91].  DSMC  simulation 
results  likely  assisted  in  decisions  regarding  the  descent  trajectory  and  orientation.  Figure  23  shows 
streamlines  around  Mir  from  SMILE  calculations  for  a  simulated  altitude  of  1 10  km. 
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Figure  24  shows  the  surface  pressure  distribution  on  Mir  due  to  Space  Shuttle  Reaction  Control  System 
(RCS)  exhaust  gases  [40],  as  part  of  simulations  using  the  DAC  DSMC  code  [43]  in  support  of  a  planned 
Space  Shuttle  docking.  DSMC  simulations  provided  detailed  information  on  plume  impingement  effects 
associated  with  various  RCS  thruster  combinations. 

Simulations  of  high  altitude  and  exo-atmospheric  rocket  exhaust  plumes  have  often  involved  an  uncoupled 
two-step  approach.  In  this  approach,  a  Navier-Stokes  CFD  code  is  employed  efficient  simulation  of 
continuum  nearfield  plume  regions  (possibly  including  portions  of  the  nozzle),  and  a  continuum 
breakdown  surface  within  the  Navier-Stokes  domain  is  then  used  as  an  inflow  boundary  for  a  subsequent 
DSMC  simulation.  DSMC  inflow  boundary  conditions  are  based  on  Navier-Stokes  simulation  results 
along  the  breakdown  surface,  and  an  accurate  solution  is  assumed  if  the  bulk  flow  is  uniformly  supersonic 
in  the  surface-normal  direction  and  directed  into  the  DSMC  domain.  This  approach  was  used  to  generate 
the  pressure  distribution  in  Figure  24.  Figure  25  shows  centerline  density  contours  from  another  example 
of  this  approach,  also  involving  the  DAC  code  [43],  for  which  interacting  plumes  are  simulated  from 
simultaneous  firings  of  three  Space  Shuttle  RCS  thrusters  [40].  Black  lines  in  Figure  25  indicate  inflow 
boundaries  used  in  the  DSMC  simulation,  and  nearfield  plume  results  within  the  regions  enclosed  by  these 
lines  are  from  an  uncoupled  Navier-Stokes  simulation.  Note  the  complex  farfield  plume  interactions 
shown  in  the  figure,  with  a  high-density  mixing  layer  bounded  on  both  sides  by  oblique  shocks.  The 
placement  of  DSMC  inflow  boundaries  along  continuum  breakdown  surfaces  in  this  type  of  flow  problem 
is  further  illustrated  in  Figure  26.  This  figure  shows  a  breakdown  surface  in  green,  based  on  calculation  of 
Bird’s  expansion  flow  parameter  [4]  from  Navier-Stokes  simulation  results,  for  a  generic  rocket  exhaust 
plume  at  high  altitude  [92].  The  Blue  surface  in  Figure  26  is  a  smoothed  approximation  of  the  continuum 
breakdown  surface,  and  is  employed  as  a  nonuniform  inflow  boundary  for  a  DSMC  simulation. 

Figure  27  shows  contours  of  surface  pressure  from  a  DSMC  simulation  of  plume  impingement  on  the 
Hubble  Space  Telescope,  due  to  Space  Shuttle  airlock  venting  during  a  planned  servicing  mission  [40]. 
Following  unanticipated  effects  of  airlock  venting  during  a  previous  Hubble  servicing  mission,  extensive 
DSMC-based  analysis  was  performed  as  part  of  mission  preparation.  The  Space  Shuttle  payload  bay, 
airlock  and  cargo,  shown  in  Figure  27,  were  also  included  in  the  simulations. 

Figure  28,  also  taken  from  DSMC  results  in  Ref.  [40],  shows  the  surface  heat  flux  distribution  over  the 
Mars  Global  Surveyor  during  a  planned  aerobraking  maneuver  in  the  upper  atmosphere  of  Mars  [93].  This 
mission  constituted  the  first  use  of  aerobraking  as  the  primary  means  of  orbit  adjustment,  and  DSMC 
simulations  played  an  important  role  in  mission  planning.  As  discussed  in  Ref.  [93],  DSMC  was  also  used 
extensively  for  post-launch  aerothermal  analysis,  following  an  unplanned  change  in  aerobraking 
configuration  due  to  an  unlatched  solar  panel. 

4.0  SUMMARY 

A  general  overview  of  rarefied  effects  in  hypersonic  flows  has  been  presented,  with  emphasis  on 
numerical  modeling,  relevant  applications  and  critical  research  issues.  Phenomena  associated  with  rarefied 
hypersonic  and  high  enthalpy  gas  flows  was  described,  including  both  translational  and  internal  energy 
nonequilibrium.  An  overview  was  provided  for  several  different  techniques  used  to  simulate  rarefied  gas 
flows,  and  representative  examples  were  discussed  for  various  techniques.  The  direct  simulation  Monte 
Carlo  (DSMC)  method,  which  is  the  most  commonly  used  simulation  method  for  high  Mach  number 
rarefied  flows,  has  been  described  in  particular  detail.  A  number  of  DSMC  examples  were  presented, 
including  detailed  comparisons  with  experimental  data  and  three  dimensional  simulations  of  complex 
engineering  problems  related  to  the  Space  Shuttle,  Mir  and  other  spacecraft  and  atmospheric  entry 
vehicles.  For  DSMC,  as  for  discrete  velocity  Boltzmann  and  other  methods,  various  shortcomings,  recent 
developments  and  areas  of  current  research  were  discussed. 
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While  any  brief  discussion  of  such  a  broad  range  of  topics  is  necessarily  incomplete,  it  is  hoped  that  the 
material  presented  here  can  serve  as  a  general  introduction  to  a  range  of  techniques  and  issues  associated 
with  simulation  of  nonequilibrium  hypersonic  flows.  It  should  be  emphasized  that  material  discussed  here 
represents  only  a  small  sample  of  the  extensive  work  in  rarefied  gas  simulation  over  the  past  several 
decades,  and  the  references  provide  far  more  detail  on  much  of  this  work. 
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Figure  1 :  Validity  range  of  different  transport  equations. 
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T  =  6000  K 


Figure  2:  Population  weighted  transition  rates  of  reactive  and  non-reactive  energy  exchanges  in 
N2  dissociation  at  6000  K;  (a)  Rate  set  1  [18, 19],  (b)  Rate  set  2  (SSH  rates)  [20]. 


Figure  3:  Variation  of  single  and  multiquantum  rates  for 
N2(i=25)+N2  ->  N2  (i’=1 5,23,24)+  N2  [21,  22,  23]. 
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Figure  4:  Translational  temperature  behind  shock  wave  for  single  and  multiquantum  VT  rates, 

Mach  19.82  nitrogen  flow,  27  Pa,  300  K. 


Figure  5:  Population  distributions  behind  shock  wave  for  single  and  multiquantum  VT  rates, 

Mach  19.82  nitrogen  flow,  27  Pa,  300  K. 
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Figure  6:  Velocity  distribution  function  inside  a  Mach  5  shock  wave  with 
hard  sphere  collision  model,  (a)  x/X  =  -10,  (b)  x/X  =  0,  (c)  x/X  =  0.8,  (d)  x/X  =  3. 


Figure  7:  Velocity  grid  resolution:  Probability  distribution  function  in  a  Mach  5  shock  wave  with 
hard  sphere  collision  model  on  coarse  and  fine  grids,  (a)  x/X  =  -1,  x/X  =  -3,  (c)  x/X  =  5. 
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Figure  8:  Accuracy  of  gas  kinetic  scheme:  Macroscopic  parameters  inside  a  Mach  3  shock  wave 
with  hard  sphere  collision  model,  (a)  density  and  temperature,  (b)  viscous  stress  xx-component, 

(c)  heat  flux  x-component. 

- Gas  Kinetic  Scheme  - Boltzmann  [33] 


Figure  9:  Accuracy  of  gas  kinetic  scheme:  Macroscopic  parameters  inside  a  Mach  5  shock  wave 
with  hard  sphere  collision  model,  (a)  density,  (b)  temperature,  (c)  streamwise  and  normal 
components  of  temperature. - Gas  Kinetic  Scheme  - Boltzmann  [33] 


1  -30 


RTO-EN-AVT-194 


NATO 

_ 1 _ 

OTAN 

Review  of  Rarefied  Gas  Effects  in  Hypersonic  Applications 


Figure  10:  Time  evolution  of  translational  temperature  components  (or  measures  of  velocity 
component  variances)  (left  figure)  and  velocity  component  covariances  (right  figure)  for  a  single 
species  undergoing  homogeneous  translational  relaxation.  Comparison  shows  good  agreement 
between  quadrature-based  approach  and  DSMC. 
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Figure  11 :  Comparison  of  density  profiles  predicted  by  collisional  Boltzmann  equation  and  BGK 
model  in  a  1 D  shock  tube,  at  t  =  0  (left)  and  t  =  0.5  (right). 
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Figure  12:  Flowchart  for  DSMC  simulation  procedures  during  each  time  step. 


Figure  13:  Variation  in  the  relative  population  of  CO  at  vibrational  quantum  level  4  as  a  function 
of  transit  time  through  a  Mach  7  shock  wave.  Red  line  on  left  shows  traditional  DSMC  result;  red 
line  on  right  shows  result  using  state  resolved  energy  exchange  procedure  [62]. 
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Figure  14:  Comparison  of  normalized  density  (left)  and  surface  heat  flux  coefficient  (right) 
between  MONACO  DSMC  and  experiment  for  hypersonic  flows  over  a  blunted  70°  cone  at  0° 
angle  of  attack  [41].  Density  is  for  Mach  20.2,  Kn  =  0.03;  heat  flux  is  for  Mach  15.6,  Kn  =  0.002. 


Figure  15:  Pressure  contours  and  streamlines  (left),  and  comparison  of  DSMC  and  experimental 
data  for  normalized  density  (right)  for  Mach  20.2  flow  of  N2  over  a  blunted  70°  cone  at  10°  angle 

of  attack,  with  Kn  =  0.03  [66]. 
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Figure  16:  Pressure  contours  and  streamlines  (left),  comparison  of  DSMC  and  experimental  heat 
flux  values  (right)  for  Mach  15.6  flow  over  double  cone,  with  Kn  =  0.002  [85]. 


Figure  17:  Translational  temperature,  surface  heat  flux  (left,  in  K  and  W/m2  respectively),  Mach 
number  and  streamlines  (right)  for  Mach  14  air  flow  over  viscous-optimized  waverider  at  100  km 

altitude  with  Kn  =  0.02. 


1  -34 


RTO-EN-AVT-194 


^■NATQ 
^  OTAN 


Review  of  Rarefied  Gas  Effects  in  Hypersonic  Applications 


p,  N/m2 


x,  m 


Figure  18:  Pressure  contours  and  streamlines  for  ballute  atmospheric  entry  configuration  [87]. 


.j 

Figure  19:  Contours  of  Mach  number  for  Mach  10  flow  of  N2  over  Apollo  command  module 
geometry  at  30°  angle  of  attack,  with  Kn  =  0.067  [88]. 
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Figure  20:  Surface  pressure  for  X-38  reentry  vehicle  at  400,000  ft  altitude  [40]. 


Figure  21 :  Density  contours  for  Space  Shuttle  Orbiter  at  170  km  altitude  [89]. 
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Figure  22:  Overall  temperature  contours  for  Beagle2  entry  into  Mars  atmosphere,  with  11°  angle 

of  attack  and  Kn  =  0.0291  [90]. 


Figure  23:  Streamlines  from  DSMC  calculation  of  flow  around  Mir  space  station  during 
controlled  descent,  at  110  km  altitude  [91]. 
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Figure  24:  Surface  pressure  on  Mir  space  station  from  Space  Shuttle  RCS  exhaust  during 

docking  maneuver  [40]. 


Figure  25:  Centerline  density  distribution  for  uncoupled  DSMC/CFD  simulations  of  Space  Shuttle 
Orbiter  during  simultaneous  firing  of  RCS  thruster  in  nose  and  two  RCS  thrusters  in  tail  [40]. 
Black  lines  indicate  DSMC/CFD  domain  boundaries. 
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Figure  26:  Continuum  breakdown  surface  (green)  and  DSMC  inflow  boundary  (blue)  used  for 
two-step  CFD/DSMC  simulation  of  rocket  exhaust  plume  at  high  altitude  [92]. 


Figure  27:  Pressure  distribution  over  Hubble  Space  Telescope  due  to  venting  from  Space 
Shuttle  airlock  during  servicing  mission  [40]. 


Figure  28:  Surface  heat  flux  over  Mars  Global  Surveyor  during  aerobraking  in  Mars  upper 

atmosphere  [40]. 
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